Analysis of simulated vs observed runoff

#Graph observed vs simulated runoff for Max Multiplicative Change Factor (MCF)

observed<- read.csv(here("climate_change_data","SWMM_results", "Mar14_observed.csv"))
simulated_max_march14 <- read.csv(here("climate_change_data","SWMM_results", "MCF_max_Mar14_simulated.csv"))
simulated_min_march14 <- read.csv(here("climate_change_data","SWMM_results", "MCF_min_Mar14_simulated.csv"))

# merge MCF max with observed
discharge_MCF_max<- merge(observed, simulated_max_march14, by = "time_step")

#merge MCF min with observed
discharge_MCF_min<- merge(observed, simulated_min_march14, by = "time_step")


summary(discharge_MCF_max)

# statical analyses
# dischargettest<- t.test(discharge_MCF_max$discharge_obs_cfs, discharge_MCF_max$simulated_flow_cfs)
# 
# regress<- lm(discharge_obs_cfs ~ simulated_flow_cfs, data = discharge_MCF_max)
# 
# regress
# 
# sim<- discharge_MCF_max$simulated_flow_cfs
# obs<- discharge_MCF_max$discharge_obs_cfs
# 
# 
# sutcliffe<- NSE(sim, obs, na.rm=TRUE, FUN=NULL, epsilon=c("Pushpalatha2012"))
# 
# sutcliffe

# write.csv(discharge,file = "discharge_obv_sim1.csv", row.names = FALSE)

graph_max<- ggplot(discharge_MCF_max, aes(x=discharge_obs_cfs, simulated_flow_cfs))+
  geom_point()

graph_max
graph_min<- ggplot(discharge_MCF_min, aes(x=discharge_obs_cfs, simulated_flow_cfs))+
  geom_point()

graph_min
calibrate_graph_max<- discharge_MCF_max %>% 
  ggplot()+
  geom_line(aes(x=time_step, y=discharge_obs_cfs), color="#000000", size = 2.5)+
  geom_line(aes(x=time_step, y=simulated_flow_cfs), color= "#009E73", size = 2.5)+
  theme_classic()+
  labs(x="Time (hours)", y="Discharge (cfs)")+
  scale_y_continuous(limits= c(0,850), breaks= seq(0,850, by= 100),expand= c(0,0))+
  scale_x_continuous(limits= c(0,12), breaks= seq(0,12, by= 2),expand= c(0,0))+
  annotate("text", label= "Observed", x=10.5, y=85, size=8)+
  annotate("text", label= "Simulated", x=10.5, y=300, size=8) +
  theme(text = element_text(size=22)) 

calibrate_graph_max  
calibrate_graph_min<- discharge_MCF_min %>% 
  ggplot()+
  geom_line(aes(x=time_step, y=discharge_obs_cfs), color="#000000", size =2.5) +
  geom_line(aes(x=time_step, y=simulated_flow_cfs), color= "#009E73", size =2.5)+
  theme_classic()+
  labs(x="Time (hours)", y="Discharge (cfs)")+
  scale_y_continuous(limits= c(0,400), breaks= seq(0,350, by= 50),expand= c(0,0))+
  scale_x_continuous(limits= c(0,12), breaks= seq(0,12, by= 2),expand= c(0,0))+
  annotate("text", label= "Simulated", x=10.5, y=45, size=8)+
  annotate("text", label= "Observed", x=10.5, y=185, size=8) +
    theme(text = element_text(size=22)) 
 

calibrate_graph_min 
# save graphs
 # ggsave('discharge_mar14_max_MCF.png', calibrate_graph_max, width = 16, height = 9, units = "in")

# ggsave('discharge_mar14_min_MCF.png', calibrate_graph_min, width = 16, height = 9, units = "in")


## combine all three lines onto one graph for easabilty 
calibrate_graph <- ggplot()+
  geom_line(data = discharge_MCF_min, aes(x=time_step, y=discharge_obs_cfs), color="#000000", size =2.5) +
  geom_line(data= discharge_MCF_min, aes(x=time_step, y=simulated_flow_cfs), color= "#009E73", size =2.5)+
  geom_line(data= discharge_MCF_max, aes(x=time_step, y=simulated_flow_cfs), color= "#D55E00", size = 2.5) +
  theme_classic()+
  labs(x="Time (hours)", y="Discharge (cfs)")+
  scale_y_continuous(limits= c(0,850), breaks= seq(0,850, by= 100),expand= c(0,0))+
  scale_x_continuous(limits= c(0,12), breaks= seq(0,12, by= 2),expand= c(0,0))+
  annotate("text", label= "Observed", x=10.5, y=200, size=8)+
  annotate("text", label= "Max MCF Simulated", x=10.5, y=655, size=8) +
   annotate("text", label= "Min MCF Simulated", x=10.5, y=55, size=8) + 
  theme(text = element_text(size=22)) 

calibrate_graph
# ggsave('discharge_mar14_MCF.png', calibrate_graph, width = 16, height = 9, units = "in")

Analysis of MCF max storm results

subcatch_max_mcf<- read_csv("subcatchments_all.csv") %>% 
  mutate(subcatchment= OBJECTID_1) %>% 
  select(subcatchment, Curve_Number, Slope, percent_imp, Area_sqft)
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   OBJECTID_1 = col_double(),
##   Area_sqft = col_double(),
##   Curve_Number = col_double(),
##   Slope = col_double(),
##   percent_imp = col_double()
## )
swm_results_max_mcf<-read_csv("Wailupe_MCF_Max_Mar14.csv") 
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   subcatchment = col_double(),
##   total_precip_in = col_double(),
##   total_runon_in = col_double(),
##   total_evap_in = col_double(),
##   total_infil_in = col_double(),
##   imperv_runoff_in = col_double(),
##   perv_runoff_in = col_double(),
##   total_runoff_in = col_double(),
##   total_runoff_10_6_gal = col_double(),
##   peak_runoff_cfs = col_double(),
##   runoff_coeff = col_double()
## )
#Characterize results by urbanization level
results_max_mcf<- merge(subcatch_max_mcf, swm_results_max_mcf, by = "subcatchment") %>% 
  mutate(runoff_normalized= 
           total_runoff_in/Area_sqft) %>% 
  mutate(Urbanization_level= 
           case_when(
             percent_imp <15 |percent_imp == 15 ~ "Natural (less than 15 % Impervious)",
             percent_imp >15 & percent_imp <45 ~ "Urbanized (Between 15 and 45 % Impervious)",
             percent_imp >44.9999 ~ "Very urbanized (More than 45 % Impervious)"
             )
           )

write.csv(results_max_mcf, file = "results_max_mcf.csv") ##Export as .csv for use with graph maps

## normalize peak discharge (cfs) for max MCF by dividing by subcatchment area 
results_max_mcf_normalized <- results_max_mcf %>% 
  mutate(peak_runoff_cfs_norm = peak_runoff_cfs/Area_sqft)

results_maxmcf_peakflow_table <- results_max_mcf %>% 
  select(subcatchment, peak_runoff_cfs, Curve_Number, Slope, percent_imp, Area_sqft) %>% 
  filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>% 
  arrange(desc(peak_runoff_cfs)) %>% 
  top_n(20,peak_runoff_cfs) %>% 
  kable(col.names=c("Subcatchment","Peak Runoff (cfs), Max MCF", 
                                "Curve Number", "Slope", 
                                "Percent Impervious", "Area (sqft)")) %>% 
  kable_styling(bootstrap_options = "striped")

results_maxmcf_peakflow_table
Subcatchment Peak Runoff (cfs), Max MCF Curve Number Slope Percent Impervious Area (sqft)
63 106.39 55.59826 11.3508475 59.595082 1260906.0
11 105.54 55.59506 0.9475857 60.176841 1461121.3
28 99.65 58.22901 9.6951530 45.648739 1400355.0
74 90.09 45.63212 0.5130342 49.951879 1593473.1
42 80.56 49.94638 12.7741916 52.075244 1097149.0
89 79.13 58.82177 12.2672167 54.379880 929708.8
78 69.55 42.51236 0.6625196 45.869531 1272102.3
20 64.99 68.12691 28.7382083 5.926749 3260351.4
62 62.88 47.11242 4.9693069 57.621143 816823.7
23 59.82 58.75298 2.0917640 60.158489 712111.1
58 57.25 46.82201 13.4737125 50.489651 794612.3
7 56.98 62.46735 24.3708501 14.101788 2174275.2
76 52.96 61.94713 12.6689327 21.946154 1533842.8
79 48.58 42.64164 0.7809627 55.674414 703208.6
75 47.57 57.16000 3.1037990 49.969260 646779.6
91 47.33 61.92057 7.2257735 49.197766 543493.1
10 38.42 59.50282 13.0349991 46.827316 527259.5
52 38.27 46.04250 14.2092259 49.068201 521285.7
67 37.74 59.21465 13.0945141 65.315445 376961.0
34 35.34 50.45073 6.5716791 55.845327 439337.8
results_maxmcf_peakflow_table_normalized <- results_max_mcf_normalized %>% 
  select(subcatchment, peak_runoff_cfs, Curve_Number, Slope, percent_imp, Area_sqft,peak_runoff_cfs_norm ) %>% 
  filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>% 
  arrange(desc(peak_runoff_cfs_norm)) %>% 
  top_n(20,peak_runoff_cfs_norm) %>% 
  kable(col.names=c("Subcatchment","Peak Runoff (cfs), Max MCF", 
                                "Curve Number", "Slope", 
                                "Percent Impervious", "Area (sqft)", "Normalized Peak Runoff (cfs/sqft)")) %>% 
  kable_styling(bootstrap_options = "striped")

results_maxmcf_peakflow_table_normalized
Subcatchment Peak Runoff (cfs), Max MCF Curve Number Slope Percent Impervious Area (sqft) Normalized Peak Runoff (cfs/sqft)
5 4.51 59.94459 4.973908 59.22529 42008.07 0.0001074
60 13.90 58.26651 10.200056 66.27512 130703.86 0.0001063
68 31.68 64.12348 16.778001 56.40240 306366.86 0.0001034
45 21.10 46.68192 10.837372 75.81188 205273.48 0.0001028
46 21.33 57.85258 9.180995 59.97949 209977.29 0.0001016
47 16.51 44.12992 9.017865 73.67874 163180.85 0.0001012
67 37.74 59.21465 13.094514 65.31544 376961.02 0.0001001
40 7.68 46.96258 7.774446 65.26444 78272.63 0.0000981
71 15.57 58.22480 1.229756 63.33945 161251.17 0.0000966
14 23.52 59.93120 10.062202 51.97367 244581.39 0.0000962
22 10.23 52.60335 5.224538 57.36459 106889.43 0.0000957
69 16.14 61.92279 23.710094 39.74707 171176.11 0.0000943
65 29.11 59.79150 12.285943 54.98557 309782.87 0.0000940
29 5.84 46.60531 4.303893 63.75090 62214.52 0.0000939
31 21.49 58.70699 10.270774 47.00970 229731.45 0.0000935
59 22.41 53.74048 9.355130 58.61615 241399.55 0.0000928
82 5.99 52.36590 12.984541 50.03241 65510.58 0.0000914
51 17.68 45.48179 3.205696 65.70618 194382.01 0.0000910
54 27.29 46.30916 11.985041 65.90850 301206.44 0.0000906
21 7.00 41.71682 2.133592 64.38962 78832.74 0.0000888
results_maxmcf_runoffcoef_table <- results_max_mcf %>% 
  select(subcatchment, runoff_coeff, Curve_Number, Slope, percent_imp, Area_sqft) %>% 
  filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>% 
  arrange(desc(runoff_coeff)) %>% 
  top_n(20,runoff_coeff) %>% 
  kable(col.names=c("Subcatchment","Runoff Coefficient, Max MCF", 
                                "Curve Number", "Slope", 
                                "Percent Impervious", "Area (sqft)")) %>% 
  kable_styling(bootstrap_options = "striped")

results_maxmcf_runoffcoef_table
Subcatchment Runoff Coefficient, Max MCF Curve Number Slope Percent Impervious Area (sqft)
60 0.821 58.26651 10.200056 66.27512 130703.86
45 0.815 46.68192 10.837372 75.81188 205273.48
68 0.813 64.12348 16.778001 56.40240 306366.86
5 0.809 59.94459 4.973908 59.22529 42008.07
67 0.809 59.21465 13.094514 65.31544 376961.02
47 0.799 44.12992 9.017865 73.67874 163180.85
71 0.793 58.22480 1.229756 63.33945 161251.17
46 0.789 57.85258 9.180995 59.97949 209977.29
65 0.770 59.79150 12.285943 54.98557 309782.87
14 0.764 59.93120 10.062202 51.97367 244581.39
40 0.762 46.96258 7.774446 65.26444 78272.63
23 0.759 58.75298 2.091764 60.15849 712111.05
91 0.755 61.92057 7.225773 49.19777 543493.09
59 0.751 53.74048 9.355130 58.61615 241399.55
54 0.748 46.30916 11.985041 65.90850 301206.44
89 0.748 58.82177 12.267217 54.37988 929708.78
29 0.747 46.60531 4.303893 63.75090 62214.52
51 0.746 45.48179 3.205696 65.70618 194382.01
63 0.746 55.59826 11.350848 59.59508 1260905.98
22 0.745 52.60335 5.224538 57.36459 106889.43
#Perform a linear regression for the SWMM max MCF storm results
results_max_mcf_regression<- lm(total_runoff_in ~Curve_Number + Slope + percent_imp + 
                              Area_sqft, data = results_max_mcf)

#Graph the relationship bw simulated runoff and impervious cover by urbanization level
runoff_imp_graph_max_mcf<- results_max_mcf %>% 
  ggplot(aes(x=percent_imp, y=total_runoff_in))+
  geom_point(aes(color=Urbanization_level))+
  labs(x= "Percent Impervious of Subcatchment", y= "Total Simulated Runoff (inches)")+
  scale_y_continuous(limits= c(0,8), breaks= seq(0,8, by= 1),expand= c(0,0.08))+
  scale_x_continuous(limits= c(0,80), breaks= seq(0,80, by= 10),expand= c(0,0))+
  scale_color_manual(name= "Urbanization Level", values= c("darkgreen", "darkseagreen", 
                                                           "darkgoldenrod1"))+
  theme_classic()

runoff_imp_graph_max_mcf

#Save runoff vs. impervious graph
# ggsave("runoff_imp_max_mcf.pdf", width = 8, height =4)
# ggsave("runoff_imp_max_mcf.png", width = 8, height =4)

#Create a table for the regression results for max mcf storm 
regress_table_max_mcf<- stargazer(results_max_mcf_regression, type ="html", digits= 2,
                          dep.var.labels = "Total Runoff (Inches)",
                       covariate.labels = c("Curve Number", "Slope", "Percent Impervious",
                                            "Area (sqft)", "Y-Intercept"),
                       omit.stat = c("rsq"))
## 
## <table style="text-align:center"><tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td><em>Dependent variable:</em></td></tr>
## <tr><td></td><td colspan="1" style="border-bottom: 1px solid black"></td></tr>
## <tr><td style="text-align:left"></td><td>Total Runoff (Inches)</td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Curve Number</td><td>0.03<sup>**</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.01)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Slope</td><td>0.05<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.01)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Percent Impervious</td><td>0.05<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.01)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Area (sqft)</td><td>0.0000<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.0000)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Y-Intercept</td><td>-0.34</td></tr>
## <tr><td style="text-align:left"></td><td>(0.61)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Observations</td><td>97</td></tr>
## <tr><td style="text-align:left">Adjusted R<sup>2</sup></td><td>0.48</td></tr>
## <tr><td style="text-align:left">Residual Std. Error</td><td>0.67 (df = 92)</td></tr>
## <tr><td style="text-align:left">F Statistic</td><td>23.33<sup>***</sup> (df = 4; 92)</td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"><em>Note:</em></td><td style="text-align:right"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>
## </table>
regress_table_max_mcf
##  [1] ""                                                                                                                                                                                                                                   
##  [2] "<table style=\"text-align:center\"><tr><td colspan=\"2\" style=\"border-bottom: 1px solid black\"></td></tr><tr><td style=\"text-align:left\"></td><td><em>Dependent variable:</em></td></tr>"                                      
##  [3] "<tr><td></td><td colspan=\"1\" style=\"border-bottom: 1px solid black\"></td></tr>"                                                                                                                                                 
##  [4] "<tr><td style=\"text-align:left\"></td><td>Total Runoff (Inches)</td></tr>"                                                                                                                                                         
##  [5] "<tr><td colspan=\"2\" style=\"border-bottom: 1px solid black\"></td></tr><tr><td style=\"text-align:left\">Curve Number</td><td>0.03<sup>**</sup></td></tr>"                                                                        
##  [6] "<tr><td style=\"text-align:left\"></td><td>(0.01)</td></tr>"                                                                                                                                                                        
##  [7] "<tr><td style=\"text-align:left\"></td><td></td></tr>"                                                                                                                                                                              
##  [8] "<tr><td style=\"text-align:left\">Slope</td><td>0.05<sup>***</sup></td></tr>"                                                                                                                                                       
##  [9] "<tr><td style=\"text-align:left\"></td><td>(0.01)</td></tr>"                                                                                                                                                                        
## [10] "<tr><td style=\"text-align:left\"></td><td></td></tr>"                                                                                                                                                                              
## [11] "<tr><td style=\"text-align:left\">Percent Impervious</td><td>0.05<sup>***</sup></td></tr>"                                                                                                                                          
## [12] "<tr><td style=\"text-align:left\"></td><td>(0.01)</td></tr>"                                                                                                                                                                        
## [13] "<tr><td style=\"text-align:left\"></td><td></td></tr>"                                                                                                                                                                              
## [14] "<tr><td style=\"text-align:left\">Area (sqft)</td><td>0.0000<sup>***</sup></td></tr>"                                                                                                                                               
## [15] "<tr><td style=\"text-align:left\"></td><td>(0.0000)</td></tr>"                                                                                                                                                                      
## [16] "<tr><td style=\"text-align:left\"></td><td></td></tr>"                                                                                                                                                                              
## [17] "<tr><td style=\"text-align:left\">Y-Intercept</td><td>-0.34</td></tr>"                                                                                                                                                              
## [18] "<tr><td style=\"text-align:left\"></td><td>(0.61)</td></tr>"                                                                                                                                                                        
## [19] "<tr><td style=\"text-align:left\"></td><td></td></tr>"                                                                                                                                                                              
## [20] "<tr><td colspan=\"2\" style=\"border-bottom: 1px solid black\"></td></tr><tr><td style=\"text-align:left\">Observations</td><td>97</td></tr>"                                                                                       
## [21] "<tr><td style=\"text-align:left\">Adjusted R<sup>2</sup></td><td>0.48</td></tr>"                                                                                                                                                    
## [22] "<tr><td style=\"text-align:left\">Residual Std. Error</td><td>0.67 (df = 92)</td></tr>"                                                                                                                                             
## [23] "<tr><td style=\"text-align:left\">F Statistic</td><td>23.33<sup>***</sup> (df = 4; 92)</td></tr>"                                                                                                                                   
## [24] "<tr><td colspan=\"2\" style=\"border-bottom: 1px solid black\"></td></tr><tr><td style=\"text-align:left\"><em>Note:</em></td><td style=\"text-align:right\"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>"
## [25] "</table>"

Analysis of MCF min storm results

subcatch_min_mcf<- read_csv("subcatchments_all.csv") %>% 
  mutate(subcatchment= OBJECTID_1) %>% 
  select(subcatchment, Curve_Number, Slope, percent_imp, Area_sqft)
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   OBJECTID_1 = col_double(),
##   Area_sqft = col_double(),
##   Curve_Number = col_double(),
##   Slope = col_double(),
##   percent_imp = col_double()
## )
swm_results_min_mcf<-read_csv("Wailupe_MCF_Min_Mar14.csv") 
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   subcatchment = col_double(),
##   total_precip_in = col_double(),
##   total_runon_in = col_double(),
##   total_evap_in = col_double(),
##   total_infil_in = col_double(),
##   imperv_runoff_in = col_double(),
##   perv_runoff_in = col_double(),
##   total_runoff_in = col_double(),
##   total_runoff_10_6_gal = col_double(),
##   peak_runoff_cfs = col_double(),
##   runoff_coeff = col_double()
## )
#Characterize results by urbanization level
results_min_mcf<- merge(subcatch_min_mcf, swm_results_min_mcf, by = "subcatchment") %>% 
  mutate(runoff_normalized= 
           total_runoff_in/Area_sqft) %>% 
  mutate(Urbanization_level= 
           case_when(
             percent_imp <15 |percent_imp == 15 ~ "Natural (less than 15 % Impervious)",
             percent_imp >15 & percent_imp <45 ~ "Urbanized (Between 15 and 45 % Impervious)",
             percent_imp >44.9999 ~ "Very urbanized (More than 45 % Impervious)"
             )
           )

write.csv(results_min_mcf, file = "results_min_mcf.csv") ##Export as .csv for use with graph maps


## normalize peak discharge (cfs) for min MCF by dividing by subcatchment area 
results_min_mcf_normalized <- results_min_mcf %>% 
  mutate(peak_runoff_cfs_norm = peak_runoff_cfs/Area_sqft)

results_minmcf_peakflow_table <- results_min_mcf %>% 
  select(subcatchment, peak_runoff_cfs, Curve_Number, Slope, percent_imp, Area_sqft) %>% 
  filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>% 
  arrange(desc(peak_runoff_cfs)) %>% 
  top_n(20,peak_runoff_cfs) %>% 
  kable(col.names=c("Subcatchment","Peak Runoff (cfs), Min MCF", 
                                "Curve Number", "Slope", 
                                "Percent Impervious", "Area (sqft)")) %>% 
  kable_styling(bootstrap_options = "striped")

results_minmcf_peakflow_table
Subcatchment Peak Runoff (cfs), Min MCF Curve Number Slope Percent Impervious Area (sqft)
63 15.29 55.59826 11.3508475 59.59508 1260906.0
28 13.56 58.22901 9.6951530 45.64874 1400355.0
42 12.15 49.94638 12.7741916 52.07524 1097149.0
11 11.91 55.59506 0.9475857 60.17684 1461121.3
89 10.85 58.82177 12.2672167 54.37988 929708.8
74 10.07 45.63212 0.5130342 49.95188 1593473.1
62 9.66 47.11242 4.9693069 57.62114 816823.7
58 8.76 46.82201 13.4737125 50.48965 794612.3
78 8.38 42.51236 0.6625196 45.86953 1272102.3
23 8.09 58.75298 2.0917640 60.15849 712111.1
76 7.40 61.94713 12.6689327 21.94615 1533842.8
7 6.79 62.46735 24.3708501 14.10179 2174275.2
75 6.65 57.16000 3.1037990 49.96926 646779.6
79 6.63 42.64164 0.7809627 55.67441 703208.6
91 5.86 61.92057 7.2257735 49.19777 543493.1
52 5.66 46.04250 14.2092259 49.06820 521285.7
10 5.46 59.50282 13.0349991 46.82732 527259.5
67 5.35 59.21465 13.0945141 65.31544 376961.0
34 5.31 50.45073 6.5716791 55.84533 439337.8
55 4.37 52.00222 14.0124563 43.23872 456026.8
results_minmcf_peakflow_normalized_table <- results_min_mcf_normalized %>% 
  select(subcatchment, peak_runoff_cfs, Curve_Number, Slope, percent_imp, Area_sqft,peak_runoff_cfs_norm ) %>% 
  filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>% 
  arrange(desc(peak_runoff_cfs_norm)) %>% 
  top_n(20,peak_runoff_cfs_norm) %>% 
  kable(col.names=c("Subcatchment","Peak Runoff (cfs), Min MCF", 
                                "Curve Number", "Slope", 
                                "Percent Impervious", "Area (sqft)", "Normalized Peak Runoff (cfs/sqft)")) %>% 
  kable_styling(bootstrap_options = "striped")

results_minmcf_peakflow_normalized_table
Subcatchment Peak Runoff (cfs), Min MCF Curve Number Slope Percent Impervious Area (sqft) Normalized Peak Runoff (cfs/sqft)
45 3.41 46.68192 10.837372 75.81188 205273.48 0.0000166
47 2.66 44.12992 9.017865 73.67874 163180.85 0.0000163
60 1.92 58.26651 10.200056 66.27512 130703.86 0.0000147
40 1.13 46.96258 7.774446 65.26444 78272.63 0.0000144
51 2.80 45.48179 3.205696 65.70618 194382.01 0.0000144
54 4.32 46.30916 11.985041 65.90850 301206.44 0.0000143
21 1.12 41.71682 2.133592 64.38962 78832.74 0.0000142
67 5.35 59.21465 13.094514 65.31544 376961.02 0.0000142
29 0.88 46.60531 4.303893 63.75090 62214.52 0.0000141
71 2.21 58.22480 1.229756 63.33945 161251.17 0.0000137
49 4.29 44.95489 11.626618 61.33458 317060.05 0.0000135
46 2.80 57.85258 9.180995 59.97949 209977.29 0.0000133
13 2.77 37.24907 2.801558 60.72206 210147.80 0.0000132
5 0.55 59.94459 4.973908 59.22529 42008.07 0.0000131
59 3.13 53.74048 9.355130 58.61615 241399.55 0.0000130
35 2.11 36.38429 2.502320 58.10591 163977.29 0.0000129
22 1.36 52.60335 5.224538 57.36459 106889.43 0.0000127
36 1.31 36.48504 3.413503 57.37835 103198.24 0.0000127
44 3.09 46.81877 9.792575 56.58471 246676.73 0.0000125
68 3.83 64.12348 16.778001 56.40240 306366.86 0.0000125
results_minmcf_runoffcoef_table <- results_min_mcf %>% 
  select(subcatchment, runoff_coeff, Curve_Number, Slope, percent_imp, Area_sqft) %>% 
  filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>% 
  arrange(desc(runoff_coeff)) %>% 
  top_n(20,runoff_coeff) %>% 
  kable(col.names=c("Subcatchment","Runoff Coefficient, Min MCF", 
                                "Curve Number", "Slope", 
                                "Percent Impervious", "Area (sqft)")) %>% 
  kable_styling(bootstrap_options = "striped")

results_minmcf_runoffcoef_table
Subcatchment Runoff Coefficient, Min MCF Curve Number Slope Percent Impervious Area (sqft)
45 0.594 46.68192 10.837372 75.81188 205273.48
47 0.586 44.12992 9.017865 73.67874 163180.85
60 0.531 58.26651 10.200056 66.27512 130703.86
40 0.528 46.96258 7.774446 65.26444 78272.63
51 0.516 45.48179 3.205696 65.70618 194382.01
54 0.514 46.30916 11.985041 65.90850 301206.44
21 0.511 41.71682 2.133592 64.38962 78832.74
29 0.510 46.60531 4.303893 63.75090 62214.52
67 0.509 59.21465 13.094514 65.31544 376961.02
71 0.493 58.22480 1.229756 63.33945 161251.17
5 0.489 59.94459 4.973908 59.22529 42008.07
49 0.484 44.95489 11.626618 61.33458 317060.05
46 0.481 57.85258 9.180995 59.97949 209977.29
13 0.473 37.24907 2.801558 60.72206 210147.80
59 0.465 53.74048 9.355130 58.61615 241399.55
22 0.463 52.60335 5.224538 57.36459 106889.43
35 0.461 36.38429 2.502320 58.10591 163977.29
36 0.458 36.48504 3.413503 57.37835 103198.24
63 0.454 55.59826 11.350848 59.59508 1260905.98
23 0.452 58.75298 2.091764 60.15849 712111.05
#Perform a linear regression for the SWMM max MCF storm results
results_min_mcf_regression<- lm(total_runoff_in ~Curve_Number + Slope + percent_imp + 
                              Area_sqft, data = results_min_mcf)

#Graph the relationship bw simulated runoff and impervious cover by urbanization level
runoff_imp_graph_min_mcf<- results_min_mcf %>% 
  ggplot(aes(x=percent_imp, y=total_runoff_in))+
  geom_point(aes(color=Urbanization_level))+
  labs(x= "Percent Impervious of Subcatchment", y= "Total Simulated Runoff (inches)")+
  scale_y_continuous(limits= c(0,7), breaks= seq(0,7, by= 1),expand= c(0,0.08))+
  scale_x_continuous(limits= c(0,80), breaks= seq(0,80, by= 10),expand= c(0,0))+
  scale_color_manual(name= "Urbanization Level", values= c("darkgreen", "darkseagreen", 
                                                           "darkgoldenrod1"))+
  theme_classic()

runoff_imp_graph_min_mcf

#Save runoff vs. impervious graph
# ggsave("runoff_imp_min_mcf.pdf", width = 8, height =4)
# ggsave("runoff_imp_min_mcf.png", width = 8, height =4)

#Create a table for the regression results
regress_table_min_mcf<- stargazer(results_min_mcf_regression, type ="html", digits= 2,
                          dep.var.labels = "Total Runoff (Inches)",
                       covariate.labels = c("Curve Number", "Slope", "Percent Impervious",
                                            "Area (sqft)", "Y-Intercept"),
                       omit.stat = c("rsq"))
## 
## <table style="text-align:center"><tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td><em>Dependent variable:</em></td></tr>
## <tr><td></td><td colspan="1" style="border-bottom: 1px solid black"></td></tr>
## <tr><td style="text-align:left"></td><td>Total Runoff (Inches)</td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Curve Number</td><td>0.001</td></tr>
## <tr><td style="text-align:left"></td><td>(0.001)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Slope</td><td>0.002<sup>*</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.001)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Percent Impervious</td><td>0.01<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.0005)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Area (sqft)</td><td>0.00<sup>**</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.00)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Y-Intercept</td><td>-0.13<sup>**</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.05)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Observations</td><td>97</td></tr>
## <tr><td style="text-align:left">Adjusted R<sup>2</sup></td><td>0.90</td></tr>
## <tr><td style="text-align:left">Residual Std. Error</td><td>0.06 (df = 92)</td></tr>
## <tr><td style="text-align:left">F Statistic</td><td>220.89<sup>***</sup> (df = 4; 92)</td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"><em>Note:</em></td><td style="text-align:right"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>
## </table>
regress_table_min_mcf
##  [1] ""                                                                                                                                                                                                                                   
##  [2] "<table style=\"text-align:center\"><tr><td colspan=\"2\" style=\"border-bottom: 1px solid black\"></td></tr><tr><td style=\"text-align:left\"></td><td><em>Dependent variable:</em></td></tr>"                                      
##  [3] "<tr><td></td><td colspan=\"1\" style=\"border-bottom: 1px solid black\"></td></tr>"                                                                                                                                                 
##  [4] "<tr><td style=\"text-align:left\"></td><td>Total Runoff (Inches)</td></tr>"                                                                                                                                                         
##  [5] "<tr><td colspan=\"2\" style=\"border-bottom: 1px solid black\"></td></tr><tr><td style=\"text-align:left\">Curve Number</td><td>0.001</td></tr>"                                                                                    
##  [6] "<tr><td style=\"text-align:left\"></td><td>(0.001)</td></tr>"                                                                                                                                                                       
##  [7] "<tr><td style=\"text-align:left\"></td><td></td></tr>"                                                                                                                                                                              
##  [8] "<tr><td style=\"text-align:left\">Slope</td><td>0.002<sup>*</sup></td></tr>"                                                                                                                                                        
##  [9] "<tr><td style=\"text-align:left\"></td><td>(0.001)</td></tr>"                                                                                                                                                                       
## [10] "<tr><td style=\"text-align:left\"></td><td></td></tr>"                                                                                                                                                                              
## [11] "<tr><td style=\"text-align:left\">Percent Impervious</td><td>0.01<sup>***</sup></td></tr>"                                                                                                                                          
## [12] "<tr><td style=\"text-align:left\"></td><td>(0.0005)</td></tr>"                                                                                                                                                                      
## [13] "<tr><td style=\"text-align:left\"></td><td></td></tr>"                                                                                                                                                                              
## [14] "<tr><td style=\"text-align:left\">Area (sqft)</td><td>0.00<sup>**</sup></td></tr>"                                                                                                                                                  
## [15] "<tr><td style=\"text-align:left\"></td><td>(0.00)</td></tr>"                                                                                                                                                                        
## [16] "<tr><td style=\"text-align:left\"></td><td></td></tr>"                                                                                                                                                                              
## [17] "<tr><td style=\"text-align:left\">Y-Intercept</td><td>-0.13<sup>**</sup></td></tr>"                                                                                                                                                 
## [18] "<tr><td style=\"text-align:left\"></td><td>(0.05)</td></tr>"                                                                                                                                                                        
## [19] "<tr><td style=\"text-align:left\"></td><td></td></tr>"                                                                                                                                                                              
## [20] "<tr><td colspan=\"2\" style=\"border-bottom: 1px solid black\"></td></tr><tr><td style=\"text-align:left\">Observations</td><td>97</td></tr>"                                                                                       
## [21] "<tr><td style=\"text-align:left\">Adjusted R<sup>2</sup></td><td>0.90</td></tr>"                                                                                                                                                    
## [22] "<tr><td style=\"text-align:left\">Residual Std. Error</td><td>0.06 (df = 92)</td></tr>"                                                                                                                                             
## [23] "<tr><td style=\"text-align:left\">F Statistic</td><td>220.89<sup>***</sup> (df = 4; 92)</td></tr>"                                                                                                                                  
## [24] "<tr><td colspan=\"2\" style=\"border-bottom: 1px solid black\"></td></tr><tr><td style=\"text-align:left\"><em>Note:</em></td><td style=\"text-align:right\"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>"
## [25] "</table>"

Maps of SWMM Results

Code setup - Load packages

library(sf) #For shapefiles
library(tmap) #For mapmaking
library(tmaptools) #For mapmaking
library(janitor) #For cleaning names

Maps for MCF max storm hotspots

results_max_mcf<- read_csv("results_max_mcf.csv") ##read in file from code above

#Combine subcatchments outline with wet storm results
subcatch_max <- read_sf(dsn = here("climate_change_data", "SWMM_results","shapefiles"), layer = "subcatch_outline") %>%
  st_transform(st_crs(4326)) %>% 
  clean_names() %>% 
  select(subcatchment = objectid_1) %>% 
  merge(results_max_mcf) %>% 
   filter(subcatchment != "5")

#  Total volume hotspots  
hotspots_max_mcf_total <- tm_basemap("OpenStreetMap.Mapnik") +
  tm_shape(subcatch_max, unit = "Miles") +
  tm_polygons("runoff_coeff", alpha = 0.8, palette = "Blues", style = "cont", n=8, 
              legend.hist = TRUE, title = "Runoff Coefficient") +
  tm_layout(title = "March 2009 Max MCF storm", inner.margins=c(.05, .05, 0.1, .53), 
            legend.position =  c(.6,.32), legend.title.size = 1.4, legend.text.size = 1) +
  tm_text("subcatchment", size = 0.3) +
  tm_scale_bar(position = c(.6,.59), breaks = c(0, 0.2, 0.4, 0.6, 0.8,1)) +
  tm_compass(position = c(.58,.52))

# tmap_save(hotspots_max_mcf_total, here("climate_change_data", "output_maps","hotspots_mcf_max.png")) ##  not exporting


#  Peak flow hotspots
hotspots_max_mcf_peak <- tm_basemap("OpenStreetMap.Mapnik") +
  tm_shape(subcatch_max, unit = "Miles") +
  tm_polygons("peak_runoff_cfs", alpha = 0.75, palette = "Greens", style = "cont", n=8, 
              legend.hist = TRUE, title = "Peak Discharge (cfs)") +
  tm_layout(title = "March 2009 Max MCF storm", inner.margins=c(.05, .05, 0.1, .53), 
            legend.position =  c(.6,.27), legend.title.size = 1.4, legend.text.size = 1) +
  tm_text("subcatchment", size = 0.3) +
  tm_scale_bar(position = c(.6,.54), breaks = c(0, 0.2, 0.4, 0.6, 0.8,1)) +
  tm_compass(position = c(.58,.47))
#tmap_save(hotspots_max_mcf_peak, here("climate_change_data", "output_maps","hotspots_max_mcf_peak.png"))

subcatch_max_peak <- subcatch_max %>% 
  filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>% 
  ggplot() +
  geom_sf(aes(fill = peak_runoff_cfs)) +
  theme_bw() + 
  scale_fill_gradient(low = "#D1F2EB", high = "#148F77", name = "Peak Discharge (cfs)",
                      breaks = c(0, 25, 50, 75, 100 )) +
  labs(title = "Max MCF") +
  geom_sf_text(aes(label = subcatchment), colour = "black", fontface = "bold", size = 1.8) +
  # geom_sf_label(aes(label = subcatchment), size = 1.5) +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
theme(axis.text.x = element_text(angle = 90))+ 
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank(),
        panel.border=element_blank()) +
  theme_void() # for faculty presentation, remove lines and axis elements 

subcatch_max_peak

# ggsave('subcatch_max_peak_removedsubs_themevoid.png', subcatch_max_peak, width = 5, height = 7, units = "in")

subcatch_max_total <- subcatch_max %>%
  filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>% 
  ggplot() +
  geom_sf(aes(fill = runoff_coeff)) +
  theme_bw() + 
  scale_fill_gradient(low = "#D6EAF8", high = "#2874A6", name = "Runoff Coefficient",
                      breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1.0),
                      limits = c(0,1.0)) +
  labs(title = "Max MCF") +
  geom_sf_text(aes(label = subcatchment), colour = "black", fontface = "bold", size = 1.8) +
  # geom_sf_label(aes(label = subcatchment), size = 1.5) +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
theme(axis.text.x = element_text(angle = 90)) +
    theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank(),
        panel.border=element_blank()) +
  theme_void()

subcatch_max_total

  #  ggsave('subcatch_max_total.png', subcatch_max_total, width = 5, height = 7, units = "in")

Maps for MCF min storm hotspots

results_min_mcf<- read_csv("results_min_mcf.csv") ##read in file from code above

#Combine subcatchments outline with wet storm results
subcatch_min <- read_sf(dsn = here("climate_change_data", "SWMM_results","shapefiles"), layer = "subcatch_outline") %>%
  st_transform(st_crs(4326)) %>% 
  clean_names() %>% 
  select(subcatchment = objectid_1) %>% 
  merge(results_min_mcf) %>% 
   filter(subcatchment != "5")

#  Total volume hotspots  
hotspots_min_mcf_total <- tm_basemap("OpenStreetMap.Mapnik") +
  tm_shape(subcatch_min, unit = "Miles") +
  tm_polygons("runoff_coeff", alpha = 0.8, palette = "Blues", style = "cont", n=8, 
              legend.hist = TRUE, title = "Runoff Coefficient") +
  tm_layout(title = "March 2009 Min MCF storm", inner.margins=c(.05, .05, 0.1, .53), 
            legend.position =  c(.6,.32), legend.title.size = 1.4, legend.text.size = 1) +
  tm_text("subcatchment", size = 0.3) +
  tm_scale_bar(position = c(.6,.59), breaks = c(0, 0.2, 0.4, 0.6, 0.8,1)) +
  tm_compass(position = c(.58,.52))
##tmap_save(hotspots_wet_total, here("5.Results_Maps", "output_maps","hotspots_wet_total.png"))


#  Peak flow hotspots
hotspots_min_mcf_peak <- tm_basemap("OpenStreetMap.Mapnik") +
  tm_shape(subcatch_min, unit = "Miles") +
  tm_polygons("peak_runoff_cfs", alpha = 0.75, palette = "Greens", style = "cont", n=8, 
              legend.hist = TRUE, title = "Peak Discharge (cfs)") +
  tm_layout(title = "March 2009 Min MCF storm", inner.margins=c(.05, .05, 0.1, .53), 
            legend.position =  c(.6,.27), legend.title.size = 1.4, legend.text.size = 1) +
  tm_text("subcatchment", size = 0.3) +
  tm_scale_bar(position = c(.6,.54), breaks = c(0, 0.2, 0.4, 0.6, 0.8,1)) +
  tm_compass(position = c(.58,.47))
##tmap_save(hotspots_wet_peak, here("5.Results_Maps", "output_maps","hotspots_wet_peak.png"))

subcatch_min_peak <- subcatch_min %>%
  filter(!subcatchment %in% c(1, 2,94,3,38, 24))  %>% 
  ggplot() +
  geom_sf(aes(fill = peak_runoff_cfs)) +
  theme_classic()+ 
  scale_fill_gradient(low = "#D1F2EB", high = "#148F77", name = "Peak Discharge (cfs)",
                      breaks = c(0, 3, 6, 9, 12, 15)) +
  labs(title = "Min MCF") +
  geom_sf_text(aes(label = subcatchment), colour = "black", fontface = "bold", size = 1.8) +
  # geom_sf_label(aes(label = subcatchment), size = 1.5) +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
theme(axis.text.x = element_text(angle = 90)) + 
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank(),
        panel.border=element_blank()) +
  theme_void() # for facily presentation, remove lines and axis elements 

subcatch_min_peak

 # ggsave('subcatch_min_peak_removedsub_themevoid.png', subcatch_min_peak, width = 5, height = 7, units = "in")

subcatch_min_total <- subcatch_min %>%
  filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>% 
  ggplot() +
  geom_sf(aes(fill = runoff_coeff)) +
  theme_bw() + 
  scale_fill_gradient(low = "#D6EAF8", high = "#2874A6", name = "Runoff Coefficient",
                      breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1.0),
                      limits = c(0,1.0)) +
  labs(title = "Min MCF") +
  geom_sf_text(aes(label = subcatchment), colour = "black", fontface = "bold", size = 1.8) +
  # geom_sf_label(aes(label = subcatchment), size = 1.5) +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
theme(axis.text.x = element_text(angle = 90)) +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank(),
        panel.border=element_blank()) +
  theme_void()

subcatch_min_total

# ggsave('subcatch_min_total.png', subcatch_min_total, width = 5, height = 7, units = "in")

Normalized Peak Discharge Maps

subcatch_max_normalized <- read_sf(dsn = here("climate_change_data", "SWMM_results","shapefiles"), layer = "subcatch_outline") %>%
  st_transform(st_crs(4326)) %>% 
  clean_names() %>% 
  select(subcatchment = objectid_1) %>% 
  merge(results_max_mcf_normalized) %>% 
   filter(subcatchment != "5")

## normalized Max MCF Peak Discharge
subcatch_max_peak_normalized <- subcatch_max_normalized %>% 
  filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>% 
  ggplot() +
  geom_sf(aes(fill = peak_runoff_cfs_norm)) +
  theme_bw() + 
  scale_fill_gradient(low = "#D1F2EB", high = "#148F77", name = "Normalized Peak Discharge (cfs/sqft)") +
  labs(title = "Max MCF") +
  geom_sf_text(aes(label = subcatchment), colour = "black", fontface = "bold", size = 1.8) +
  # geom_sf_label(aes(label = subcatchment), size = 1.5) +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
theme(axis.text.x = element_text(angle = 90))+ 
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank(),
        panel.border=element_blank()) +
  theme_void() # for faculty presentation, remove lines and axis elements 

subcatch_max_peak_normalized

# ggsave('subcatch_max_peak_removedsubs_themevoid_norm.png', subcatch_max_peak_normalized, width = 5, height = 7, units = "in")

subcatch_min_normalized <- read_sf(dsn = here("climate_change_data", "SWMM_results","shapefiles"), layer = "subcatch_outline") %>%
  st_transform(st_crs(4326)) %>% 
  clean_names() %>% 
  select(subcatchment = objectid_1) %>% 
  merge(results_min_mcf_normalized) %>% 
   filter(subcatchment != "5")

## normalized Min MCF Peak Discharge
subcatch_min_peak_normalized <- subcatch_min_normalized %>%
  filter(!subcatchment %in% c(1, 2,94,3,38, 24))  %>% 
  ggplot() +
  geom_sf(aes(fill = peak_runoff_cfs_norm)) +
  theme_classic()+ 
  scale_fill_gradient(low = "#D1F2EB", high = "#148F77", name = "Normalized Peak Discharge (cfs/sqft)") +
  labs(title = "Min MCF") +
  geom_sf_text(aes(label = subcatchment), colour = "black", fontface = "bold", size = 1.8) +
  # geom_sf_label(aes(label = subcatchment), size = 1.5) +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
theme(axis.text.x = element_text(angle = 90)) + 
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank(),
        panel.border=element_blank()) +
  theme_void() # for facily presentation, remove lines and axis elements 

subcatch_min_peak_normalized

 # ggsave('subcatch_min_peak_removedsub_themevoid_norm.png', subcatch_min_peak_normalized, width = 5, height = 7, units = "in")

binned runoff coefficient maps for min and max MCF

library(RColorBrewer)
library(metR)


# View a single RColorBrewer palette by specifying its name
display.brewer.pal(n = 7, name = 'Blues')

# Hexadecimal color specification
brewer.pal(n = 7, name = "Blues")
## [1] "#EFF3FF" "#C6DBEF" "#9ECAE1" "#6BAED6" "#4292C6" "#2171B5" "#084594"
mycols <- brewer.pal(n = 11, name = "Blues")
mybreaks_max<- c(0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9)

subcatch_min_total <- subcatch_min %>%
  filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>%
  mutate(runoff_coeff_cut = cut_number(runoff_coeff, n = 10)) %>%
  ggplot() +
  geom_sf(aes(fill = runoff_coeff_cut)) +
  theme_bw()  +
  scale_fill_gradientn(colours = mycols,
                       breaks= mybreaks_max,
    super = metR::ScaleDiscretised,
    name = "Runoff Coefficient") +
  theme(legend.position = "bottom") +
  labs(title = "Min MCF") +
  geom_sf_text(aes(label = subcatchment), colour = "black", fontface = "bold", size = 3.1) +
  # geom_sf_label(aes(label = subcatchment), size = 1.5) +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
theme(axis.text.x = element_text(angle = 90)) +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank(),
        panel.border=element_blank()) +
  theme_void()

subcatch_min_total

  # ggsave('subcatch_min_total_binned.png', subcatch_min_total, width = 5, height = 7, units = "in")

mycols <- brewer.pal(n = 11, name = "Blues")
mybreaks_max<- c(0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9)

subcatch_max_total <- subcatch_max %>%
  filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>%
  mutate(runoff_coeff_cut = cut_number(runoff_coeff, n = 10)) %>%
  ggplot() +
  geom_sf(aes(fill = runoff_coeff_cut)) +
  theme_bw()  +
  scale_fill_gradientn(colours = mycols,
                       breaks= mybreaks_max,
    super = metR::ScaleDiscretised,
    name = "Runoff Coefficient") +
  labs(title = "Max MCF") +
  geom_sf_text(aes(label = subcatchment), colour = "black", fontface = "bold", size = 3.1) +
  # geom_sf_label(aes(label = subcatchment), size = 1.5) +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
theme(axis.text.x = element_text(angle = 90)) +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank(),
        panel.border=element_blank()) +
  theme_void()

subcatch_max_total

 # ggsave('subcatch_max_total_binned.png', subcatch_max_total, width = 5, height = 7, units = "in")

binned peak flow maps for min and max MCF

mycols2 <- brewer.pal(n = 11, name = "Greens")
mybreaks2_min<- c(0.0, 2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0)

subcatch_min_peak <- subcatch_min %>%
  filter(!subcatchment %in% c(1, 2,94,3,38, 24))  %>%
  mutate(peak_runoff_cfs_cut = cut_number(peak_runoff_cfs, n = 11)) %>% 
  ggplot() +
  geom_sf(aes(fill = peak_runoff_cfs_cut)) +
  theme_classic()+
  scale_fill_gradientn(colours = mycols2,
                       breaks= mybreaks2_min,
    super = metR::ScaleDiscretised,
    name = "Peak Discharge (cfs)") +
  labs(title = "Min MCF") +
  geom_sf_text(aes(label = subcatchment), colour = "black", fontface = "bold", size = 3.1) +
  # geom_sf_label(aes(label = subcatchment), size = 1.5) +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
theme(axis.text.x = element_text(angle = 90)) + 
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank(),
        panel.border=element_blank()) +
  theme_void() # for facily presentation, remove lines and axis elements 

subcatch_min_peak

 # ggsave('subcatch_min_peak_removedsub_themevoid_binned.png', subcatch_min_peak, width = 5, height = 7, units = "in")

mycols3 <- brewer.pal(n = 11, name = "Greens")
mybreaks3_max<- c(0.0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100)

subcatch_max_peak <- subcatch_max %>% 
  filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>%
  mutate(peak_runoff_cfs_cut = cut_number(peak_runoff_cfs, n = 11)) %>% 
  ggplot() +
  geom_sf(aes(fill = peak_runoff_cfs_cut)) +
  theme_classic()+
  scale_fill_gradientn(colours = mycols3,
                       breaks= mybreaks3_max,
    super = metR::ScaleDiscretised,
    name = "Peak Discharge (cfs)") +
  labs(title = "Max MCF") +
  geom_sf_text(aes(label = subcatchment), colour = "black", fontface = "bold", size = 3.1) +
  # geom_sf_label(aes(label = subcatchment), size = 1.5) +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
theme(axis.text.x = element_text(angle = 90))+ 
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank(),
        panel.border=element_blank()) +
  theme_void() # for faculty presentation, remove lines and axis elements 

subcatch_max_peak

# ggsave('subcatch_max_peak_removedsubs_themevoid_binned.png', subcatch_max_peak, width = 5, height = 7, units = "in")

combine max and min runoff coefficient results into single table showing the catchment #, min, and max runoff coefficient, slope, percent impervious, and area (sqft)

results_min_mcf_rc <- results_min_mcf %>%
  select(subcatchment, Slope, percent_imp, Area_sqft, runoff_coeff) %>%
  filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>%
  arrange(desc(runoff_coeff)) %>%
  top_n(20,runoff_coeff)

results_max_mcf_rc <- results_max_mcf %>%
  select(subcatchment, Slope, percent_imp, Area_sqft, runoff_coeff) %>%
  filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>%
  arrange(desc(runoff_coeff)) %>%
  top_n(20,runoff_coeff)

results_runoffcoef_table <- results_min_mcf_rc %>%
  full_join(results_max_mcf_rc, by=c("subcatchment", "Slope", "percent_imp", "Area_sqft")) %>%
  kable(col.names=c("Subcatchment","Slope", "Percent Impervious", "Area (sqft)", "Runoff Coefficient, Min MCF",
                                "Runoff Coefficient, Max MCF")) %>%
  kable_styling(bootstrap_options = "striped")

results_runoffcoef_table
Subcatchment Slope Percent Impervious Area (sqft) Runoff Coefficient, Min MCF Runoff Coefficient, Max MCF
45 10.837372 75.81188 205273.48 0.594 0.815
47 9.017865 73.67874 163180.85 0.586 0.799
60 10.200056 66.27512 130703.86 0.531 0.821
40 7.774446 65.26444 78272.63 0.528 0.762
51 3.205696 65.70618 194382.01 0.516 0.746
54 11.985041 65.90850 301206.44 0.514 0.748
21 2.133592 64.38962 78832.74 0.511 NA
29 4.303893 63.75090 62214.52 0.510 0.747
67 13.094514 65.31544 376961.02 0.509 0.809
71 1.229756 63.33945 161251.17 0.493 0.793
5 4.973908 59.22529 42008.07 0.489 0.809
49 11.626618 61.33458 317060.05 0.484 NA
46 9.180995 59.97949 209977.29 0.481 0.789
13 2.801558 60.72206 210147.80 0.473 NA
59 9.355130 58.61615 241399.55 0.465 0.751
22 5.224538 57.36459 106889.43 0.463 0.745
35 2.502320 58.10591 163977.29 0.461 NA
36 3.413503 57.37835 103198.24 0.458 NA
63 11.350848 59.59508 1260905.98 0.454 0.746
23 2.091764 60.15849 712111.05 0.452 0.759
68 16.778001 56.40240 306366.86 NA 0.813
65 12.285943 54.98557 309782.87 NA 0.770
14 10.062202 51.97367 244581.39 NA 0.764
91 7.225773 49.19777 543493.09 NA 0.755
89 12.267217 54.37988 929708.78 NA 0.748

combine max and min peak runoff (cfs) results into single table showing the catchment #, min, and max peak runoff (cfs)

results_min_mcf_peak <- results_min_mcf %>%
  select(subcatchment, Slope, percent_imp, Area_sqft, peak_runoff_cfs) %>%
  filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>%
  arrange(desc(peak_runoff_cfs)) %>%
  top_n(20,peak_runoff_cfs)

results_max_mcf_peak <- results_max_mcf %>%
  select(subcatchment, Slope, percent_imp, Area_sqft, peak_runoff_cfs) %>%
  filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>%
  arrange(desc(peak_runoff_cfs)) %>%
  top_n(20,peak_runoff_cfs)

results_peakflow_table <- results_min_mcf_peak %>%
  full_join(results_max_mcf_peak, by=c("subcatchment", "Slope", "percent_imp", "Area_sqft")) %>%
  kable(col.names=c("Subcatchment","Slope", "Percent Impervious", "Area (sqft)","Peak Runoff (cfs), Min MCF",
                                "Peak Runoff (cfs), Max MCF")) %>%
  kable_styling(bootstrap_options = "striped")

results_peakflow_table
Subcatchment Slope Percent Impervious Area (sqft) Peak Runoff (cfs), Min MCF Peak Runoff (cfs), Max MCF
63 11.3508475 59.595082 1260906.0 15.29 106.39
28 9.6951530 45.648739 1400355.0 13.56 99.65
42 12.7741916 52.075244 1097149.0 12.15 80.56
11 0.9475857 60.176841 1461121.3 11.91 105.54
89 12.2672167 54.379880 929708.8 10.85 79.13
74 0.5130342 49.951879 1593473.1 10.07 90.09
62 4.9693069 57.621143 816823.7 9.66 62.88
58 13.4737125 50.489651 794612.3 8.76 57.25
78 0.6625196 45.869531 1272102.3 8.38 69.55
23 2.0917640 60.158489 712111.1 8.09 59.82
76 12.6689327 21.946154 1533842.8 7.40 52.96
7 24.3708501 14.101788 2174275.2 6.79 56.98
75 3.1037990 49.969260 646779.6 6.65 47.57
79 0.7809627 55.674414 703208.6 6.63 48.58
91 7.2257735 49.197766 543493.1 5.86 47.33
52 14.2092259 49.068201 521285.7 5.66 38.27
10 13.0349991 46.827316 527259.5 5.46 38.42
67 13.0945141 65.315445 376961.0 5.35 37.74
34 6.5716791 55.845327 439337.8 5.31 35.34
55 14.0124563 43.238717 456026.8 4.37 NA
20 28.7382083 5.926749 3260351.4 NA 64.99

comparing historical vs Max MCF runoff coeffiecients from Dornan et al., 2020

## runoff coefficients differences
dornan_results_top20_rc <- read_csv("~/Desktop/Aloha-Aina-Master/climate_change_data/SWMM_results/dornan_results_top20_rc.csv") %>% 
  select(subcatchment, runoff_coeff) # read in dorana et al RC results

results_max_mcf_rc_summary <- results_max_mcf_rc %>%  select(subcatchment, runoff_coeff) # select just our subcatchment and runoff coefficients to be able to subtract 

combined_rc_results <- dornan_results_top20_rc %>%  right_join(results_max_mcf_rc_summary, by = c("subcatchment")) # combine dornan results with our max MCF reslts 

combined_rc_results = combined_rc_results %>%  mutate(subtract = runoff_coeff.y-runoff_coeff.x) # add column with difference between max MCF rc results - dornan historical rc results with March 14, 09 storm 

## peak flow differences 
dornan_results_top20_peakflow <- read_csv("~/Desktop/Aloha-Aina-Master/climate_change_data/SWMM_results/dornan_results_top20_peakflow.csv") %>% 
  select(subcatchment, peak_runoff_cfs) # read in Team Kahuawais peak flow results

results_max_mcf_peak_summary <- results_max_mcf_peak %>%  select(subcatchment, peak_runoff_cfs) # select just our subcatchment and peak runoff cfs to be able to subtract 

combined_peak_results <- dornan_results_top20_peakflow %>%  right_join(results_max_mcf_peak_summary, by = c("subcatchment"))

combined_peak_results = combined_peak_results %>%  mutate(subtract = peak_runoff_cfs.y-peak_runoff_cfs.x)